# Replication code for 
# Weidmann, Nils B. 2011. “Violence ”from above” or ”from below”? The Role of Ethnicity in Bosnia’s Civil War.” Journal of Politics 73(4).

library(stats)
library(maptools)
library(spdep)
library(flexmix)

 # change this to wherever the file "bosnia_replication.shp" is located
setwd("/Users/nilsw/Projects/Bosnia/R/")

# read dataset
bosnia <- readShapePoly("bosnia_replication.shp",proj4string=CRS("+proj=longlat +datum=WGS84"))

# neighborhood matrix
bosnia.nb <- poly2nb(bosnia)
wmat <- nb2mat(bosnia.nb,style="W")

# generate variables
bosnia$incomepc <- 0.1* bosnia$TOTINCOME / bosnia$TOTPOP91
bosnia$roaddensity <- 1000* bosnia$roadlength/bosnia$area
bosnia$l.distserbia <- log10(bosnia$distserbia)
bosnia$l.distcroatia <- log10(bosnia$distcroati)
bosnia$borderdist <- pmin(bosnia$l.distserbia, bosnia$l.distcroatia)

# test for spatial correlation
moran.test((bosnia$KILLED/bosnia$TOTPOP91), nb2listw(bosnia.nb))


# spatial regression models (Table 1)

m1 <- lagsarlm(KILLED/TOTPOP91 ~  pol_local + borderdist + roaddensity + urban91 + incomepc, data=bosnia, listw=nb2listw(bosnia.nb))
summary(m1)

m2 <- lagsarlm(KILLED/TOTPOP91 ~  strat91 + borderdist + roaddensity + urban91 + incomepc, data=bosnia, listw=nb2listw(bosnia.nb))
summary(m2)

m3 <- lagsarlm(KILLED/TOTPOP91 ~  pol_local + strat91  + borderdist  + roaddensity + urban91 + incomepc, data=bosnia, listw=nb2listw(bosnia.nb))
summary(m3)


# Supplementary models in the online appendix

m1.a <- lagsarlm((KILLED+MISSING)/TOTPOP91 ~  pol_local + borderdist + roaddensity + urban91 + incomepc, data=bosnia, listw=nb2listw(bosnia.nb))
summary(m1.a)

m2.a <- lagsarlm((KILLED+MISSING)/TOTPOP91 ~  strat91 + borderdist + roaddensity + urban91 + incomepc, data=bosnia, listw=nb2listw(bosnia.nb))
summary(m2.a)

m3.a <- lagsarlm((KILLED+MISSING)/TOTPOP91 ~  pol_local + strat91 +  borderdist + roaddensity + urban91 + incomepc, data=bosnia, listw=nb2listw(bosnia.nb))
summary(m3.a)

#  Mixture model in Table 2

bosnia.copy <- bosnia[!is.na(bosnia$pol_local),]
model <- FLXMRglmfix(nested=list(k=c(1,1), formula=c(~ strat91, ~ pol_local)))
result <- stepFlexmix(KILLED/TOTPOP91 ~  1, k=2, model=model, concomitant=FLXPmultinom(~ l.distserbia), data=bosnia.copy, nrep=20)
summary(refit(result))
summary(refit(result), which="concomitant")

#  Mixture model in the online appendix

bosnia.copy <- bosnia[!is.na(bosnia$pol_local),]
model <- FLXMRglmfix(nested=list(k=c(1,1), formula=c(~ strat91, ~ pol_local)))
result <- stepFlexmix((KILLED+MISSING)/TOTPOP91 ~  1, k=2, model=model, concomitant=FLXPmultinom(~ l.distserbia), data=bosnia.copy, nrep=20)
summary(refit(result))
summary(refit(result), which="concomitant")





